// 版权归Go作者所有。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

import "math/bits"

func zero(x uint64) uint64 {
	if x == 0 {
		return 1
	}
	return 0
	// 无分支：
	// 返回（（x>>1 | x&1）-1）>>63 
}

func nonzero(x uint64) uint64 {
	if x != 0 {
		return 1
	}
	return 0
	// 无分支：
	// 返回1-（（x>>1 | x&1）-1）>>63 
}

func shl(u1, u2 uint64, n uint) (r1, r2 uint64) {
	r1 = u1<<n | u2>>(64-n) | u2<<(n-64)
	r2 = u2 << n
	return
}

func shr(u1, u2 uint64, n uint) (r1, r2 uint64) {
	r2 = u2>>n | u1<<(64-n) | u1>>(n-64)
	r1 = u1 >> n
	return
}

// SHRCPRESS将两个字
// 值的底部n+1位压缩为一位。结果等于向右移位n的值
// 除了结果的第0位是
// 设置为底部n+1位的按位OR之外。
func shrcompress(u1, u2 uint64, n uint) (r1, r2 uint64) {
	// TODO:这里的性能对这些分支的顺序/位置非常敏感。n==0是常见的
	// 足以进入快速路径。也许需要进行更多的测量以找到最佳的顺序/位置？
	switch {
	case n == 0:
		return u1, u2
	case n == 64:
		return 0, u1 | nonzero(u2)
	case n >= 128:
		return 0, nonzero(u1 | u2)
	case n < 64:
		r1, r2 = shr(u1, u2, n)
		r2 |= nonzero(u2 & (1<<n - 1))
	case n < 128:
		r1, r2 = shr(u1, u2, n)
		r2 |= nonzero(u1&(1<<(n-64)-1) | u2)
	}
	return
}

func lz(u1, u2 uint64) (l int32) {
	l = int32(bits.LeadingZeros64(u1))
	if l == 64 {
		l += int32(bits.LeadingZeros64(u2))
	}
	return l
}

// 拆分将b拆分为符号、有偏指数和尾数。
// 它将正常值的尾数加上隐式1位，
// 并对次正常值进行规格化。
func split(b uint64) (sign uint32, exp int32, mantissa uint64) {
	sign = uint32(b >> 63)
	exp = int32(b>>52) & mask
	mantissa = b & fracMask

	if exp == 0 {
		// 如果值低于正常值，则将其正常化。
		shift := uint(bits.LeadingZeros64(mantissa) - 11)
		mantissa <<= shift
		exp = 1 - int32(shift)
	} else {
		// 添加隐式1位
		mantissa |= 1 << 52
	}
	return
}

// FMA返回x*y+z，仅用一个舍入计算。
// （即，FMA返回x、y和z的融合乘加）
func FMA(x, y, z float64) float64 {
	bx, by, bz := Float64bits(x), Float64bits(y), Float64bits(z)

	// Inf或NaN或零。最多会出现一个舍入。
	if x == 0.0 || y == 0.0 || z == 0.0 || bx&uvinf == uvinf || by&uvinf == uvinf {
		return x*y + z
	}
	// 分别处理非有限z。计算x*y+z，其中
	// x和y是有限的，但z是无限的，应该总是得到z。
	if bz&uvinf == uvinf {
		return z
	}

	// 输入（次）正常。
	// 将x、y、z拆分为符号、指数、尾数。
	xs, xe, xm := split(bx)
	ys, ye, ym := split(by)
	zs, ze, zm := split(bz)

	// 计算乘积p=x*y作为符号、指数、两个单词尾数。
	// 从指数开始。“正常”位尚未减去。
	pe := xe + ye - bias + 1

	// pm1:pm2是产品p的双字尾数。
	// 左移以保留产品中的顶部位。实际上
	// 将106位乘积向左移位21。
	pm1, pm2 := bits.Mul64(xm<<10, ym<<11)
	zm1, zm2 := zm<<10, uint64(0)
	ps := xs ^ ys // 乘积符号

	// 规格化为62位
	is62zero := uint((^pm1 >> 62) & 1)
	pm1, pm2 = shl(pm1, pm2, is62zero)
	pe -= int32(is62zero)

	// 交换加法操作数so | p |>=| z | 
	if pe < ze || pe == ze && pm1 < zm1 {
		ps, pe, pm1, pm2, zs, ze, zm1, zm2 = zs, ze, zm1, zm2, ps, pe, pm1, pm2
	}

	// 对齐有效位
	zm1, zm2 = shrcompress(zm1, zm2, uint(pe-ze))

	// 计算结果有效位，必要时进行规格化。
	var m, c uint64
	if ps == zs {
		// 添加（pm1:pm2）+（zm1:zm2）
		pm2, c = bits.Add64(pm2, zm2, 0)
		pm1, _ = bits.Add64(pm1, zm1, c)
		pe -= int32(^pm1 >> 63)
		pm1, m = shrcompress(pm1, pm2, uint(64+pm1>>63))
	} else {
		// 减去（pm1:pm2）-（zm1:zm2）
		// TODO:我们应该取消特殊情况吗？
		pm2, c = bits.Sub64(pm2, zm2, 0)
		pm1, _ = bits.Sub64(pm1, zm1, c)
		nz := lz(pm1, pm2)
		pe -= nz
		m, pm2 = shl(pm1, pm2, uint(nz-1))
		m |= nonzero(pm2)
	}

	// 舍入并断开与偶数的关系
	if pe > 1022+bias || pe == 1022+bias && (m+1<<9)>>63 == 1 {
		// 舍入值溢出指数范围
		return Float64frombits(uint64(ps)<<63 | uvinf)
	}
	if pe < 0 {
		n := uint(-pe)
		m = m>>n | nonzero(m&(1<<n-1))
		pe = 0
	}
	m = ((m + 1<<9) >> 10) & ^zero((m&(1<<10-1))^1<<9)
	pe &= -int32(nonzero(m))
	return Float64frombits(uint64(ps)<<63 + uint64(pe)<<52 + m)
}
